Wetland to Upland Gradient
Anthony Stewart
Questions:
Does a model derived gradient between wetlands and uplands correspond to a gradient of soil development along a catena/hillslope characterized in a hydropedological framework?
If the ends of the wetland/upland classification are accurate, does the intermediate range correspond with hypothesized soil formation characteristics according to a hydropedological framework?
hydropedology hillslope model from Bailey et
al., 2014
The Hubbard Brook Experimental Forest
From Hubbard Brook website https://hubbardbrook.org/about-the-forest/:
- “The Hubbard Brook Experimental Forest is located on unceded Abenaki land in the southwestern part of the White Mountain National Forest in New Hampshire. The experimental forest occupies a 3,037-hectare (~8,700-acre) bowl-shaped valley that is oriented east–west, with a series of small watersheds occupying the valley’s north- and south-facing slopes.”
Hydropedological units in the Hubbard Brook Experimental Forest include:
- Drier more upland soils: Bhs, Typical (Podzol)
- Wetter soils: Bh, Histosol
- Intermediate/mesic soils: Bimodal, Inceptisol
## [1] "ATTRIBUTE" "WETLAND_TY" "ACRES" "Shape_Leng" "Shape_Area"
Merge hydrography and NWI
Get training data - wetland and upland points
Make terrain metrics
hbslp_3 <- SlpAsp(hbdem, w = c(3, 3), metrics = "slope", filename = "UplandWetlandGradient/data/derived_data/hbslp_3.tif",
overwrite = T)
hbslp_27 <- SlpAsp(hbdem, w = c(27, 27), metrics = "slope", filename = "UplandWetlandGradient/data/derived_data/hbslp_27.tif",
overwrite = T)
hbslp_81 <- SlpAsp(hbdem, w = c(81, 81), metrics = "slope", filename = "UplandWetlandGradient/data/derived_data/hbslp_81.tif",
overwrite = T)
hbtpi_3 <- TPI(hbdem, w = c(3, 3), shape = "rectangle", stand = "none",
na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbtpi_3.tif",
overwrite = T)
hbtpi_27 <- TPI(hbdem, w = c(27, 27), shape = "rectangle", stand = "none",
na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbtpi_27.tif",
overwrite = T)
hbtpi_81 <- TPI(hbdem, w = c(81, 81), shape = "rectangle", stand = "none",
na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbtpi_81.tif",
overwrite = T)
hbcurv_3 <- Qfit(hbdem, w = c(3, 3), metrics = c("meanc", "profc",
"planc"), unit = "degrees", na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbcurv_3.tif",
overwrite = T)
hbcurv_27 <- Qfit(hbdem, w = c(27, 27), metrics = c("meanc",
"profc", "planc"), unit = "degrees", na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbcurv_27.tif",
overwrite = T)
hbcurv_81 <- Qfit(hbdem, w = c(81, 81), metrics = c("meanc",
"profc", "planc"), unit = "degrees", na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbcurv_81.tif",
overwrite = T)Slope at different scales
TPI at different scales
Curvature at different scales
Slope at 3m vs. 81m
TPI at 3m vs. 81m
Mean Curvature at 3m vs. 81m
Is resampling the same as setting the window size larger?
## [1] "class" "dem" "slp_3" "slp_27" "slp_81"
## [6] "meancurv_3" "prof_curv_3" "plan_curv_3" "meancurv_27" "prof_curv_27"
## [11] "plan_curv_27" "meancurv_81" "prof_curv_81" "plan_curv_81" "tpi_3"
## [16] "tpi_27" "tpi_81" "x" "y" "meanc"
## [21] "profc" "planc"
## class : SpatVector
## geometry : points
## dimensions : 2246, 22 (geometries, attributes)
## extent : 274043.4, 283958.1, 4866042, 4871958 (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / UTM zone 19N (EPSG:26919)
## names : class dem slp_3 slp_27 slp_81 meancurv_3 prof_curv_3
## type : <fact> <num> <num> <num> <num> <num> <num>
## values : UPL 333.2 14.31 12.14 9.286 -0.02777 -0.04499
## UPL 745.8 15.09 8.697 12.81 -0.03448 0.00313
## UPL 584 14.21 15.18 15.11 0.001648 0.007331
## plan_curv_3 meancurv_27 prof_curv_27 (and 12 more)
## <num> <num> <num>
## -0.01056 0.0003643 -0.0001417
## -0.07208 -0.0001351 0.008627
## -0.004036 0.002487 0.003541
Set up random forest model
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
## The following object is masked from 'package:dplyr':
##
## combine
## The following object is masked from 'package:spatialEco':
##
## combine## Confusion Matrix and Statistics
##
## Reference
## Prediction UPL WET
## UPL 1264 0
## WET 0 309
##
## Accuracy : 1
## 95% CI : (0.9977, 1)
## No Information Rate : 0.8036
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 1
##
## Mcnemar's Test P-Value : NA
##
## Sensitivity : 1.0000
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 1.0000
## Prevalence : 0.8036
## Detection Rate : 0.8036
## Detection Prevalence : 0.8036
## Balanced Accuracy : 1.0000
##
## 'Positive' Class : UPL
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction UPL WET
## UPL 528 25
## WET 14 108
##
## Accuracy : 0.9422
## 95% CI : (0.9219, 0.9586)
## No Information Rate : 0.803
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.8115
##
## Mcnemar's Test P-Value : 0.1093
##
## Sensitivity : 0.9742
## Specificity : 0.8120
## Pos Pred Value : 0.9548
## Neg Pred Value : 0.8852
## Prevalence : 0.8030
## Detection Rate : 0.7822
## Detection Prevalence : 0.8193
## Balanced Accuracy : 0.8931
##
## 'Positive' Class : UPL
## 
AUC/ROC
## [1] "Area under the ROC curve"
## [[1]]
## [1] 0.8930999
Predict maps
PDF/kernel density plot
## Warning: [hist] a sample of17% of the cells was used (of which 5% was NA)
## Warning: [hist] a sample of17% of the cells was used (of which 5% was NA)
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:randomForest':
##
## margin
## # A tibble: 9 × 2
## hpu result
## <fct> <dbl>
## 1 Bh 0.148
## 2 Bhs 0.0447
## 3 Bhsm 0.0616
## 4 Bi 0.0750
## 5 E 0.0643
## 6 H 0.503
## 7 I 0.193
## 8 O 0.0280
## 9 T 0.0612
## Warning: Removed 4 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Warning: Removed 4 rows containing missing values (`geom_point()`).
Is there a model?
##
## Call:
## lm(formula = WET ~ hpu, data = hb_pedons_ext)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19318 -0.06195 -0.04470 0.00085 0.74792
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.14808 0.01844 8.029 1.55e-14 ***
## hpuBhs -0.10339 0.02941 -3.515 0.000498 ***
## hpuBhsm -0.08648 0.05249 -1.648 0.100330
## hpuBi -0.07308 0.02750 -2.657 0.008246 **
## hpuE -0.08379 0.03514 -2.385 0.017634 *
## hpuH 0.35492 0.11142 3.185 0.001577 **
## hpuI 0.04509 0.04196 1.075 0.283294
## hpuO -0.12008 0.03934 -3.052 0.002445 **
## hpuT -0.08693 0.02392 -3.634 0.000322 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1554 on 346 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.1119, Adjusted R-squared: 0.09132
## F-statistic: 5.447 on 8 and 346 DF, p-value: 1.759e-06
## Analysis of Variance Table
##
## Response: WET
## Df Sum Sq Mean Sq F value Pr(>F)
## hpu 8 1.0524 0.13155 5.4471 1.759e-06 ***
## Residuals 346 8.3558 0.02415
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Warning: Removed 1 rows containing non-finite values (`stat_ydensity()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 30 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 30 rows containing missing values (`geom_point()`).
Is there a model in W3?
##
## Call:
## lm(formula = WET ~ hpu, data = w3_pedons_ext)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.12298 -0.04011 -0.01673 0.01687 0.25955
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.116994 0.034248 3.416 0.00184 **
## hpuBhs podzol -0.066107 0.048434 -1.365 0.18244
## hpuBimodal podzol -0.087680 0.076582 -1.145 0.26130
## hpuE podzol -0.078464 0.050134 -1.565 0.12805
## hpuLithic Histosol -0.088856 0.052315 -1.698 0.09977 .
## hpuTypical podzol 0.007971 0.055224 0.144 0.88620
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09687 on 30 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1668, Adjusted R-squared: 0.02797
## F-statistic: 1.201 on 5 and 30 DF, p-value: 0.3323
## Analysis of Variance Table
##
## Response: WET
## Df Sum Sq Mean Sq F value Pr(>F)
## hpu 5 0.05637 0.0112739 1.2015 0.3323
## Residuals 30 0.28151 0.0093836##
## Call:
## lm(formula = WET ~ water_table, data = w3_pedons_ext)
##
## Residuals:
## 3 6 11 17 19 28 35
## -0.11031 -0.06748 -0.04230 0.02109 0.25246 -0.03985 -0.01360
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1561297 0.1427138 1.094 0.324
## water_table -0.0008764 0.0024205 -0.362 0.732
##
## Residual standard error: 0.13 on 5 degrees of freedom
## (30 observations deleted due to missingness)
## Multiple R-squared: 0.02555, Adjusted R-squared: -0.1693
## F-statistic: 0.1311 on 1 and 5 DF, p-value: 0.7321
## Analysis of Variance Table
##
## Response: WET
## Df Sum Sq Mean Sq F value Pr(>F)
## water_table 1 0.002215 0.0022146 0.1311 0.7321
## Residuals 5 0.084466 0.0168931